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In this report, we investigate dynamical robustness of a complex network to noise injected through one of its 
nodes. We focus on synchronization of coupled nonlinear systems and, as a special instance, we address the 
classical consensus protocol for linear integrators. We establish an exact closed-form expression of the 
synchronization error for the consensus protocol and an approximate result for chaotic units. While 
structural robustness is known to be significantly affected by attacks targeted to network hubs, our results 
posit that dynamical robustness is controlled by both the topology of the network and the dynamics of the 
units. We provide examples where hubs perform better or worse than isolated nodes. 

Communication networks, transportation infrastructures, and power grids are all subjected to failures that 
often involve a significant portion of the system despite the local nature of the initial fault. Thus, error and 
attack tolerance, cascading failures, and, in general, robustness of complex networks have been the body of 
intense research^"^^. 

To quantify network robustness, generally defined as the ability of the network to withstand accidental events, 
different measures have been proposed: conditional connectivity^ restricted connectivity^, super connectivity^, 
fault diameter^, expansion parameter^, isoperimetric number^, and natural connectivity^. Other measures of 
network robustness have been developed in the theoretical framework of statistical physics and percolation 
theory^^"^^. These efforts seek to elucidate the effect of the removal of a fraction of nodes (or links) on character- 
istic properties of a network, such as its diameter, largest component, and efficiency. 

Such investigations have demonstrated that topology is a determinant of robustness, whereby heterogeneous 
(scale-free) networks are highly robust against random attacks, which can instead severely impact homogenous 
networks. Nevertheless, targeted attacks to hubs of scale-free networks can dramatically affect network properties. 
Indeed, the collapse of the entire network can be caused by a single node whose failure propagates to neighboring 
nodes promoting cascading damages and avalanches^^"^^. This phenomenon is well explained by modeling the 
dynamics of flow redistribution on the network to unveil the mechanisms underlying the large breakdowns 
observed in real systems, such as the Internet or electrical power grids^^'^^. 

Most of the existing studies deal with the robustness of the network structure and are based on the assumption 
that the failure of a node (or a fraction of nodes) is equivalent to a complete loss of its (their) functionalities. Here, 
we take a different approach to understand the effect of a partial malfunctioning on the dynamics of the network. 
This scenario is expected to be relevant for biological networks^^, distributed sensors^^ and, in general, complex 
systems composed of coupled dynamical units^°. Notably, robustness with respect to the system dynamics has 
been recently investigated in^\ where networks of diffusively coupled second-order periodic oscillators are 
considered and node failure is modelled as the inactivation of oscillations (quenching). 

In this report, we study networks of identical coupled dynamical units (including chaotic oscillators) and 
elucidate on the effect of partial malfunctioning of a network node on synchronization. We model the phenom- 
enon of node failure by injecting noise into the dynamics of a node. In particular, we consider a white Gaussian 
signal with zero-mean and variance cr^. Increasing the noise level, that is, cr^, the synchronization error increases 
and the units may eventually leave the basin of attraction for sufficiently strong noise. Our primary objective is to 
establish a mathematical framework for assessing dynamical robustness and investigating the complex interplay 
between network topology and node dynamics. We derive tractable expression for the synchronization error and 
introduce a novel network parameter, which allows to rank the nodes in terms of their impact on the network 
dynamical robustness. We analyze four representative dynamics for coupled dynamical systems and explore 
different network topologies to illustrate the proposed methodology and demonstrate that both the topology and 
the dynamics are determinants of node ranking. 
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Results 

Evaluation of dynamical robustness to noise injection into a node. 

Here, we establish a mathematical framework for assessing dynami- 
cal robustness and investigating the interplay between network 
topology and node dynamics. In particular, we consider a network 
of N dynamical units, which synchronize in the absence of noise. 
Dynamical robustness is quantified through the synchronization 
error as a function of the noise variance. To compute such error, 
we first derive the (stochastic) dynamical equations of the transverse 
modes of the system. From these equations, we relate the synchro- 
nization error to the statistical properties of the transverse modes 
and, more in detail, to their joint moments. Based on the derived 
expression of the error, we then introduce a novel network 
parameter, which allows to rank the nodes in terms of their impact 
on the network dynamical robustness. 

We consider a network of N dynamical units, described by the 
following coupled equations: 

N 

X, = F(xO - Kj2gij^- i^j) + ^/H,f/ (1) 

;=i 

for / = 1, where x/gR"^ represents the state vector of the ith 

node, F the individual dynamics, k the coupling coefficient, the ijth 
entry of the network Laplacian G, and Hx the coupling function. The 
term rj is a. zero -mean Gaussian white noise of variance cr^ and 
HrjEU^ is a vector of zeros and ones that models how the noise 
influences the dynamics of a node. We hypothesize that the noise 
acts on a single node, say node /z, so that lifi = h and c^/ = 0 if / 

To study the effect of the noise intensity and injection site on the 
network synchronization, we suppose that the network synchronizes 
in the absence of the noise and evaluate the synchronization error 
with respect to the noise variance. We consider the following syn- 
chronization error: 



(2) 



where {')t indicates averaging with respect to Mn a defined obser- 
vation window. For a sufficiently large observation window and 
assuming the system to be ergodic, the error (see Methods) is rewrit- 
ten as: 



N 



(3) 



where E[ • ] is the expected value, T contains the right eigenvectors of 
the graph Laplacian G, and y;^ for A; = 2, AT are the transverse 
modes (see Methods). 

Now, we distinguish between a special case of synchronization, 
namely the classical consensus problem^^, for which an exact formula 
of E [yIyi] is derived, and the more general case of synchronization of 
chaotic units, for which an approximation of E [y^yj is established. 

We first consider a consensus problem, that is, X/gR and F = 0. In 
this case (see Methods), E[)/;t)//] is given by: 



^\ykyi] = 



(4) 



As expected, in (3) is linear with respect to cr and both the network 
topology and the node dynamics are important in shaping the net- 
work dynamical robustness. 

We now consider the case of synchronization of chaotic units. To 
this aim, we rewrite the dynamical equations of the transverse modes 
(see Methods) in terms of the scaled Gaussian noise rj^ = 

(Thk — N~^^jTjk)r] (with zero mean and variance {(tIY = 
(Thk — N~^ ^jTjk) (7^), and we numerically compute the evolution 



of the modes for each value of cr^ and calculate E [jJyi] as a function 
of a. 

Consistently with the consensus problem, numerical simulations 
show that E [jJyi] is proportional to cr^cr^*. Thus, we define the func- 



tion M so that M(jj^, Yi) 



Since the transverse modes 



only differ for the term = K:y^, in the following we use the notation 
M{(Xk, oci). Notably, for the consensus problem, such function 
reduces to l/(a^ + a/). In general, M depends only on the node 
dynamics, the coupling function, and the modality for which noise 
is injected into the node. 

To quantify the effect of noise injected in the network at the hth 
node, we define the normalized quantity = ^, which takes the 
following general form: 



Ph'- 



\ 



N N ( ^ 1 \ / ^1 



T^7 KM (ait, a/) 
(5) 



Injection of noise in correspondence of nodes with small p\i leads to 
smaller synchronization error, compared to the case when noise is 
applied to nodes with large p/j. Hence, the network is more robust to 
noise injected at sites with low p^. We remark that p^ depends on the 
topology of the network and on the dynamics of the unit. In contrast 
with structural robustness, which depends only on topology, we find 
that dynamics is an essential element of dynamical robustness. 

Interplay of topology and dynamics. We illustrate the concept of 
dynamical robustness through some examples. We consider two 
types of networks, that is, scale-free (SF) and Erdos-Renyi (ER) 
networks, and four different dynamics for the network units: a 
simple integrator, a Chua's circuit^^, a Rossler's system^^ and a 
Chen's system^^ (see Methods). These systems are selected on the 
basis of an exhaustive analysis on a multitude of dynamical systems 
(not reported here) to highlight the salient features of dynamical 
robustness that are unveiled by our study. We first investigate the 
agreement between our theoretical arguments and numerical or 
experimental data and then analyse py^ for some networks with a 
large number of nodes. 

We consider a SF and an ER network with AT = 100 nodes and 
average degree equal to 2 (see^^ and Methods). The coupling coef- 
ficient K is chosen so that the networks synchronize in the absence of 
noise and Ky2 does not vary. We integrate equation set (1) and 
evaluate the synchronization error (2). As expected, we find that 
E}^ is linear in a and we define = ^ to compare the theoretical 
results based on p/j with numerical findings based on ^/j. 

In Fig. 1, we report py^ vs by^ for both types of networks. The 
agreement between simulations and theoretical expectations is 
almost perfect; theoretical expectations correctly reproduce the rank- 
ing of nodes based on the network performance in the presence of 
noise. In fact, nodes can be equivalently ranked according to either 
the value of p/j or of by^. 

Next, we compare theoretical expectations with a network of real 
electronic oscillators. The experimental setup is that described in^^ 
and consists of a network of six Chua's circuits coupled according to 
the network topology shown in Fig. 2(a). Fig. 2(b) shows the com- 
parison between py^ (theoretical expectation) and (experimental 
data) for each network node, demonstrating a very good agreement 
between theoretical and experimental results. 

To offer some statistical relevance to the proposed method, we 
numerically analyze p^ for networks with a large number of nodes. In 
particular, we consider SF and ER networks with AT = 1000 nodes 
and average degree equal to 2. Fig. 3 illustrates the results for four case 
studies. In panels (a,b,c,d), we display the value of p/j ordered in 
ascending order for a SF and an ER network. In panels (e,f,g,h), we 
display the degree of node h, labelled as /c/j, vs. p/j. The network 
dynamics are those of Fig. 1. SF networks demonstrate worse 
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Figure 1 | ph vs. for each of the iV nodes of a ER network (a,b,c,d) and a SF network (e,f,g,h) with N = 100 and average degree equal to 2. The 

dynamics of a unit is: (a)-(e) integrator; (b)-(f) Chua's circuit; (c)-(g) Rossler's system; and (d)-(h) Chen's system. The coupUng coefficient is /c = 1 (a,e); 
K = 18.38 (b); K = 2.97 (c); k = 45.93 (d); k = 13.81 (f); k = 1.73 (g); and k = 34.54 (h). 



performance than ER networks for most of the nodes in the case of 
consensus and Chua's dynamics. Instead, in the case of Rossler's 
systems, the SF network performs better than the ER network with 
the exception of a small percentage of nodes, which are the network 
hubs. In the case of Chen's systems, the ER network performs better 
than the SF network for about half of the nodes. 

Based on the premise that the precise ranking of the nodes in terms 
of dynamical robustness is dictated by Ph, we observe some general 
trends in the node ranking with respect to topological measures, such 
as the node degree. Specifically, for the consensus problem and for 
the network of Chua's circuits, we find that nodes with high degree 
have better robustness to noise as compared to nodes with low 
degree. On the contrary, for networks of Rossler's systems, the 
opposite holds, whereby nodes with low degree display better robust- 
ness to noise than nodes with high degree. In the case of Chen's 
systems, the correlation between kh and ph is weaker. In fact, the 
nodes that display poor dynamical robustness are those with high 
degree. Yet, there are also several low degree nodes with limited 
performance as shown in panel (h). 

Our analysis demonstrates that robustness to noise is the result of 
the interplay between dynamics and topology, whereby dynamical 
robustness with respect to noise injected in select node cannot be 
inferred only on the basis of its degree or other topological measures 
such as, for instance, eigenvector centrality. To further elucidate this 




Figure 2 | Experimental results: (a) topology of the network of 6 Chua's 
circuits and (b) comparison between /?h (theoretical expectation) and 
(experimental data) for each node of the network. 



interplay, we report in Fig. 4 the function M{ak, oci) for the different 
systems investigated. We distinguish between functions M that are 
monotonically decreasing with respect to each of their variables 
and a/ and those that are not. Specifically, we define as 5PT^ the 
family of functions that are monotonically decreasing with respect to 
each of the two variables. Thus, for consensus and Chua's circuit 
Me^^, while for Rossler's and Chen's dynamics M^fM^ We 
observe that this difference is crucial to determine the opposite beha- 
viors found in the correlation between node degree and node robust- 
ness, for instance, for the Chua's circuit and the Rossler's system. 

The interplay between topology and dynamics in dynamical 
robustness can be further illustrated through the analysis of an exem- 
plary network, which allows for the derivation of closed-form expres- 
sions. In particular, we focus on a star network topology for which we 
can find Ph for every network node. In this network, one node, 
labelled as node 1 and referred to as the hub, is connected to all 
the other nodes of the network, referred to as the leaves. All the leaves 
have degree equal to one, while the hub has degree N — 1. The star 
network Laplacian has eigenvalues equal to: y i = 0, 72,. . 1 ^ 1 > and 

= N. For such network, lengthy, but trivial, calculations lead to the 
following formulas for pi and P2,...,n- 



These expressions are particularly simple since they depend only on 
two values of the function that is, A^(a2, 0C2) and M{(Xn, a^v). It 
is easy to note that pi is greater than P2,...,n (that is, the leaves 
perform better than the hubs) if and only if M{(Xn,ocn)> 
M{oL2, 0C2). Given the form of M for the Chen's system, 
M{(Xn, ocn) can be either greater or lower than M{(X2, (X2), which 
implies that according to the values of 1X2 and (and, ultimately, on 
the coupling strength k), the leaves can perform better than the hub 
or vice versa. This is illustrated in Fig. 5 for a star network ofN = 5 
Chen's systems. For k = Ki = 30, the hub of the network has superior 
dynamical robustness performance than peripheral nodes. In fact, 
K:iy2 and Kiy^are such that A1(K:iy^, Kiy^) <M{Kiy2, Kiy2) (panel 
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Figure 3 | Numerical results for networks with N = 1000 nodes: (a,b,c,d) and (e,f,g,h) degree of node h (ic/,). The dynamics of a unit is: (a) -(d) 
consensus; (b)-(e) Chua's circuit; (c)-(f) Rossler's system; and (d)-(h) Chen's system. The coupling coefficient k: is 1 for consensus; 14.45 and 16.22 for SF 
and ER networks of Chua's circuits; 1.81 and 2.03 for SF and ER networks of Rossler's systems; and 14.88 and 17.49 for SF and ER networks of Chen's 
systems. 



(a), black circles). Instead, for k = K2 = 60 the opposite holds, 
M{k2Yn, K2Yn)^ ^{^2Y2^ ^2^2) (panel (a), red circles), and the 
leaves outperform the hub. 

The analysis of the star network confirms that, depending on the 
shape of M, the same topology may display opposite correlations 
between node degree and node robustness. In the general case 
(shown in Fig. 1(h)), for which M ^^)Jl^, a clear correlation betwe- 
en node degree and node robustness cannot be readily established. 

Discussion 

Here, we have introduced a mathematical framework to understand 
the dynamical robustness of complex networks to noise injected into 
their nodes. The approach lends itself to the formulation of a novel 
performance metric, see Eq. (5), which allows for the quantification 
of dynamical robustness as a function of network topology and node 



dynamics. We have considered a spectrum of dynamics and topol- 
ogies to demonstrate the complex response of networks to noise 
injection. In contrast with structural robustness, where attacks tar- 
geted to hubs are generally more severe than removal of low-degree 
nodes, injection of noise into hubs may be more or less dangerous 
than injection of noise into low- degree nodes depending on the unit 
dynamics. Experimental evidence of this second scenario has been 
reported for small networks (limited to up to six nodes) in^^. A 
similar role of low-degree nodes has been also observed in^\ where 
periodic oscillators have been considered and failure has been mod- 
elled through quenching rather than injected noise. Finally, we men- 
tion that another measure of network robustness in terms of the 
dynamical properties of the network has been introduced in^°, where 
dynamic vulnerability, that is, the collective response of the network 
to a finite time perturbation into a single node, is studied for a 
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1 : 0.04 




Figure 4 | The function M.{<Xk^ oc/) for (a) consensus; (b) Chua's circuit; (c) Rossler's system; and (d) Chen's system. To better visualize the specific 
features of each function, the range for which they are reported is not the same. 
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Figure 5 | (a) The function M{(Xk, a/) for the Chen's system, (b-c) Sh vs. ^for a star network of AT = 5 nodes for two different values of the coupUng a*. In 
(b) (/r = 30), the hub of the network has better dynamical robustness performance than peripheral nodes. In (c) (a* = 60), the opposite holds. The two 

cases correspond to a different location of the eigenvalues (K:y2,3,4 = k and Ky^ = S/c), leading to two different sets of values along the curve M{oCk, oci) 
marked in (a) as black (/c = 30) or red (k = 60) circles. Note that /z = 1 is the hub and /z = 2, 5 are the leaves. 



network of Rossler's systems. The conclusion that the hubs are not 
the more vulnerable nodes is in agreement with our results. 

We have found that networks of Chua's circuits and Rossler's 
systems are representative of two opposite behaviors that can be 
exhibited by a network in terms of its dynamical robustness to noise. 
Specifically, we have demonstrated that, in networks of Chua's cir- 
cuits, the ranking of the nodes mirrors the node degree, while, in 
networks of Rossler's systems, low-degree nodes are more robust to 
noise injection. Furthermore, we have observed (as in the case of 
networks of Chen's systems) that there are scenarios in which a clear 
correlation between node degree and node robustness cannot be 
established. The different behaviors exhibited by these systems are 
explained by the shape of the function in Eq. (5). For systems with 
Mefyjl'J^, which we classify as dynamical robustness class I systems, 
we have observed a degree-correlated behavior (that is, higher degree 
nodes have higher dynamical robustness performance than lower 
degree nodes). On the other hand, for systems with M^'OXl^, which 
we classify as dynamical robustness class 11 systems, a clear correlation 
between node degree and node robustness cannot be established. For 
dynamical robustness class II systems, different scenarios may arise 
depending on the system (topological and dynamical) parameters, 
including degree-correlated and anti-degree correlated behaviors. 
Thus, our main conclusion is that dynamical robustness is affected 
by both the topology of the network and the dynamics of the units 
and that the parameter p/j is a suitable measure for node dynamical 
robustness. 

Methods 

Dynamical equations of transverse modes. Using the Kronecker algebra, Eqs. (1) 
can be rewritten in the compact form: 

x = F(x)-k;G(x)Hx(x)+S(x)H^^ (7) 

where x = [xi, X2, . . ., x^l^ H,(x) = [H,(xi), H,(x2, . . ., H,(xn)]^, and S = [c^i, (^2, • • 

Due to the presence of the noise, the oscillators cannot synchronize, that is, the 
solution Xi = . . . = XTsf is not feasible. Thus, we define the synchronization error with 
respect to a virtually synchronous state s whose evolution proxies the individual node 
dynamics when the noise is uniformly distributed in the network: s = F(s) + iV~ ^ H,^ ^7. 
We now linearize Eq. (7) around the synchronous state s to obtain 



C = [In (8) DF - k;G (X) DHx] C + S (X) ?7 



(8) 



where C = x - In® s, IjveK^is the vector of all ones, Jacobians are evaluated on the 
synchronous state, and S = [c^i — ^ , (^2 ~ ^ ? • • • ? '^n ~ ^] -To study Eq. (8), we 
define a new set of state variables y = T"^ (x) \J^, where T contains the right 
eigenvectors of the graph Laplacian G and I^ is the identity matrix of order m. For the 



sake of simplicity, we suppose that the network is undirected so that G is symmetric 
and T-^ = r. Since ^^t= [Thi-EjTji/N, ...,ThN- EjTjn/nY, we have: 



y, = [DF-^y,DHjy,+ n,-N-' ^tAh^vj 



(9) 



where for A; = 1, . . ., Nare the ordered eigenvalues of G. For k= 1, yi = 0 and Tji = 
1 for j = 1, N so that the variation along the synchronization manifold evolves 
independently of the injected noise being controlled by the Lyapunov exponents of 
the individual dynamics. For k = 2, ...,N, Eqs. (9) represent the dynamical equations 
of the transverse modes. The noise injected into the hth node affects each transverse 
mode Yk and is weighted by the term (T^^ — J2jTjk / ■ 

The synchronization error (2) is now rewritten in terms of the transverse modes. In 
fact, since ^= TiT. t (where x = iV-^ ^ /X,- is the 

average trajectory that can be approximated as x = s + y J, and = x, — s = J2k TikYk^ 
the error (2) can be rewritten as: 



Sh-- 



{\\Ti,2Y2^Ti^3Y3^ ■ ■ ■ ^Ti^NjNW^) ^ 



(10) 



Now, we illustrate the form of the dynamical equations of the transverse modes and 
the expression of E [yJyi] for the two cases investigated (the classical consensus 
problem and the more general case of synchronization of chaotic units). 
For the consensus problem (x/gIR and F = 0), Eq. (9) becomes 



(11) 



(12) 



By noting that for a continuous-time system in the form: 

ei = -ai8i-\-r]^ 

82 = —a282+r\2 

with ai, a2> 0 and rji and rj2 zero-mean Gaussian processes with 

variance Gi and (T2, respectively, E[ei£2] = ^iykyi] is given by 

ai +a2 

{Thk- Y^j-TjkjlThi- J2j^Tji 

^\ykyi] = r^^H 

In the case of synchronization of chaotic units, Eq. (9) is rewritten in terms of the 
scaled Gaussian noise rjl = (T^^ — AT"^ J2j'^jk) ^ (with zero mean and variance 
{(JlY = {Thk-N-' ZjTjk) (t'), that is, 

y, = [DF-^y,DHjy, + H,,7* (13) 

In this case, E [yJyi] is calculated from Eq. (13). In particular, we numerically 
integrate Eq. (13) with the Euler-Maruyama integration method, we compute the 
evolution of the modes for each value of al and calculate E [yJyi] ^s a function of a. 

Model equations. The network for the consensus problem is described by Eqs. (1) 
with m = 1, XiER, F{Xi) = 0, H^{Xi) = Xi and = I. 

The network of Chua's circuits is described by Eqs. (1) with m = 3, x^- = j, x, 2> 
Xi^^V, H^ixi) = 1, 0, 0]^, = [1, 0, 0]^, and F(x,) given by 
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F(x.-) = 



acfe - (mi + l^i -0.5(mo -mi)(|x,-,i + 1| - - 1|)) 



(14) 



The parameters are chosen so that the Chua's circuit is chaotic in the double scroll 
regime: Oic = 9, jSc = 14.3, thq = —1/7, and nii = 2/7. 

The network of Rossler's systems is described by Eqs. (1) with m = 3, = %f,2> 
Xi^^V, H,(x,) = [0, %,-2, 0]^ = [1, 0, 0]^ and F(x,) given by 

-10(%.-,2+^i,3) 

F(xO= 10(%,-,i+a«%,-,2) (15) 

The parameters are chosen so that the Rossler's system is chaotic: aR = 0.2, Br = 0.2, 
and Cr = 9. The Rossler's equations are scaled by a factor equal to 10 to make the time 
scale similar to that of the Chua's circuit. 

The network of Chen's systems is described by Eqs. (1) with m = 3, x, = [x, i, %, 2> 
Xi,^r, H,(x,) = [0, Xi^2, OV, = [1, 0, 0]^ and F(x,) given by 

0.5{acXi^2-acXui) 

F(xj)= 0.5((cc-ac-x/,3)%i,i +ccX/,2) (16) 

0.5{x^Xi^2-bcXu3)) 

The following parameters have been chosen: ac = 35, be = 8/3, and Cc = 28. The 
Chen's equations are scaled by a factor equal to 0.5 to render the time scale similar to 
that of the other systems investigated. 

Eqs. (1) are integrated with the Euler-Maruyama method for a window of duration 
equal to 2000 with a time step of 0.001 (which is a time step adequate for all the three 
dynamics investigated) and noise standard deviation a in the range [0, 0.5] (with a 
step of 0.02). The synchronization error (2) is evaluated on the last 500 time units of 
the simulation. 

The Master Stability Function. In our analysis, the coupling coefficient k is chosen so 
that the network synchronizes in the absence of the noise. This selection is informed 
by the Master Stability Function (MSF)^^ 

According to this approach, a block diagonalized variational equation of the form 

= [DV — Kyj^DH]^!^ represents the dynamics of the system around the synchron- 
ization manifold. Here, is the h-th eigenvalue ofG,h= 1, . . ., AT, and DF and DH are 
the Jacobian matrices of F and H computed at the synchronous state. Therefore, the 
blocks of the diagonalized variational equation differ from each other only for the 
term Ky^. To investigate the synchronization properties with respect to different 
topologies, the variational equation is studied as a function of a generic eigenvalue a 
(for the sake of simplicity, we limit the discussion to the case of undirected networks, 
for which the Laplacian has real eigenvalues). This leads to the definition of the 
Master Stability Equation (MSE): C= [DF-(xDH]C. 

The maximum (conditional) Lyapunov exponent Imax of the MSE is studied as a 
function of a, thus obtaining the MSF, that is, Q(a). Then, the stability of the syn- 
chronization manifold in a given network is evaluated by computing the eigenvalues 
yfc (with /c = 2, . . . , AT) of the matrix G and studying the sign of Q at the points = Ky^. 
If for all eigenmodes with A; = 2, . . ., N, Q is negative, then the synchronous state is 
stable at the given coupling strength k. 

In the literature, systems with always positive Q are referred to as type I MSF 
systems ^ systems having Q negative for a larger than a given constant are referred to 
as type II MSF systems^ and systems having Q negative only in a compact interval, 
that is, Q < 0 in [ai, 1X2], are referred to as type III MSF systems^ 

All the systems investigated in this report have type II MSF so that our analysis is 
not biased by the MSF type. 

Networks' construction. We consider Erdos-Renyi (ER) networks, that is, networks 
having a connectivity distribution (the probability that a given node has a given 
number of links) following a Poisson-like distribution, and scale-free (SF) networks, 
that is, networks with a heterogeneous degree distribution (a power-law distribution). 
Both networks are generated by using the model reported in^^. The model has a single 
tunable parameter which allows to interpolate between ER and scale- free SF networks 
as far as the degree distribution is concerned. 
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